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Abstract. Analyses of recent helioseismic data have produced ample evidence for substantial dynamical variation 
of the differential rotation within the solar convection zone. Given the inevitable difficulties in resolving the 
precise nature of variations at deeper layers, much effort has recently gone into determining theoretically the 
expected modes of behaviour, using nonlinear dynamo models. Two important limitations of these models are 
that they have so far included only one form of nonlinearity, and as yet they have not taken into account the 
density stratification in the solar convection zone. Here we address both of these issues by studying the effects 
of including density stratification, as well as including an a-quenching nonlinearity in addition to the previously 
studied effects of the Lorentz force on the differential rotation. We find that observationally important features 
found in the earlier uniform density models remain qualitatively unchanged, although there are quantitative 
differences. This is important as it provides more realistic theoretical predictions to be compared with and guide 
observations, especially in the deeper regions where the uncertainties in the inversions are larger. However the 
presence of an effective alpha-quenching nonlinearity significantly reduces the amplitudes of the oscillations. 

Key words. Sun: magnetic fields - torsional oscillations - activity - density stratification 



1. Introduction 

Recent analyses of the helioseismic data have produced 
a wealth of information concerning the dynamical modes 
of behaviour associated with the differential rotation in 
the solar convection zone (CZ). Previous analyses of data 
from the Michelson Doppler Imager (MDI) instrument 
on board the SOHO spacecraft (Howe et al. 2000a) and 
the Global Oscillation Network Group (GONG) project 
(Antia & Basu 2000) have found that torsional oscilla- 
tions penetrate into the convection zone, to depths of at 
least 10 percent in radius, but did not agree on the be- 
haviour in the poorly resolved region near the bottom of 
the convection zone. In order to understand these observed 
dynamical phenomena theoretically, and partly motivated 
by these apparent discrepancies, we have previously made 
detailed studies of the observed variations in the differ- 
ential rotation in the CZ in the framework of nonlinear 
dynamo models which included the nonlinear action of 
the azimuthal component of the Lorentz force of the dy- 
namo generated magnetic field on the solar angular veloc- 
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ity (Covas et al. 2000a, b; Covas et al. 2001a, b; Tavakol et 
al. 2002). These studies produced a number of important 
robust features, that may be treated as their predictions. 
These include: (i) the robust existence of torsional oscilla- 
tions that penetrated all the way down to the bottom of 
the convection zone for near-critical and moderately su- 
percritical dynamo regimes; (ii) the existence of parameter 
ranges for which the supercritical models show spatiotem- 
poral fragmentation, with different oscillatory modes of 
behaviour at the top and the bottom of the CZ; (iii) the 
existence of equatorial as well as polar branches in the an- 
gular velocity residuals; and (iv) the presence of secondary 
modulations on these equatorial branches. Predictions (i) 
and (ii) are in principle able to account for different as- 
pects of earlier observations. Interestingly, the most recent 
inversions by Vorontsov et al. (2002a) demonstrate further 
important agreement with our predictions. In particular 
they provide strong evidence that these oscillations do in 
fact penetrate down to the bottom of the CZ. Furthermore 
they find clear evidence for the presence of polar branches 
as well as of modulations on the equatorial branch (see 
Vorontsov et al. 2002b for preliminary results, and also 
Vorontsov et al. 2003). 

Despite these successes, there are two important is- 
sues concerning these models that require further study. 
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Firstly, they ignore an important feature of the real Sun, 
namely the presence of a strong density stratification in 
the solar convection zone and, secondly, they involve only 
one form of nonlinearity - that due to the Lorentz force 
acting on the large-scale velocity field (here just the ro- 
tational flow). A more realistic model would be expected 
to include other forms of nonlinearity, such as the effects 
of the Lorentz force on the small-scale motions, for ex- 
ample as naively modelled by a nonlinear a-quenching. 
Investigating the consequences of these shortcomings is 
particularly important, given the limitations on the length 
of the helioseismic data sets and the resulting uncertain- 
ties in the inversions, especially regarding the behaviour 
in the lower parts of the CZ (see e.g. Vorontsov et al 2002). 

The aim of this paper is to improve our previous model 
by (i) including a density stratification in our previous 
model, and (ii) by adding a second form of nonlinearity 
in the form of an explicit a-quenching. In contrast, in our 
previous work we imposed a radial variation of alpha that 
we felt might be consistent with a strong quenching in the 
overshoot region where the magnetic field was expected 
to be strong (by putting a = in this region). We can 
now make a detailed study of both the effects of density 
stratification and the competing forms of nonlinearity on 
various observationally important features of our model. 
Throughout the underlying zero order angular velocity is 
taken to be consistent with the recent inversions of the 
helioseismic (MDI) data. 

The structure of the paper is as follows. In the next 
section we outline our model, and Sects. and 0] contain 
our detailed results. In Sect. we attempt to draw some 
conclusions. 

2. The model 

2.1. Equations 

We write the large-scale velocity field in the form u = 
v(j> — ^V?7, and put v = vo + v', vq = Qorsmd being the 
underlying (assumed known) solar rotation law. Thus we 
solve the standard mean-field dynamo equation 



dB 

— = V x (u x B + aB - 77V x B), 



(1) 



for the mean magnetic field and the corresponding equa- 
tion 



dv' (VxBjxB, 



dt 



Airpr sinf 



4> + p- 1 V.T, 



(2) 



(r the stress tensor) for the deviations from the underly- 
ing (~ mean) rotational velocity vo', as usual 77 and v are 
turbulent transport coefficients. We proceed as described 
in Moss & Brooke (2000) and in previous papers in this 
series, Nondimensionalizing in terms of the solar radius R 
and a time R 2 /rj^, where 770 is the maximum value of 77, 
and putting Q — a — ao<5, 77 — 77077, B — BqB and 

v' = f2* Rv' . This gives dynamo parameters R a = aoR/r/Q, 
R u = il*i? 2 /?7o, Pi = vq/vo, where fl* is the solar surface 



equatorial angular velocity and B — (A47rp r] Q*) 1 / 2 , 
where A is an arbitrary scaling factor and po is the density 
at the lower boundary of the model. As discussed in Moss 
& Brooke (2000), the physical field B B is independent 
of the value of A - for 'historical' reasons we set A = 10 4 . 
P r is the turbulent Prandtl number. 

The model includes a density stratification discussed 
in Section l2~2l below, as well as a second source of nonlin- 
earity in the form of an a - quenching discussed in Section 
Ol below. 

In this investigation, f2o is given in 0.64 < r < 1 by an 
interpolation on the MDI data obtained from 1996 to 1999 
(Howe et al. 2000a). For a we put a = a T (r)f(0), where 
f{9) was chosen to be sin 2 cos#. This form has been used 
previously (sec e.g. Rudigcr & Brandenburg 1995) and 
its choice here is simply to make the butterfly diagrams 
more realistic. For a r we took the following choices. In the 
absence of a-quenching, a r was set to be a r = 1 for 0.7 < 
r < 0.8 with cubic interpolation to zero at r = r$ and r = 
1. In presence of a-quenching, we took a r = 1 throughout 
the convection zone. We use the convention that a r > 
and R a < 0, to get the correct sense of migration of field 
patterns. Also, in order to take into token account the 
likely decrease in the turbulent diffusion coefficient 77 in 
the overshoot region, we allowed a simple linear decrease 
from 77 = 1 at r = 0.8 to 77 = 0.5 in r < 0.7. 

At the surface r = R we took vacuum boundary condi- 
tions, whereby the poloidal field within r — R is smoothly 
joined, by a matrix multiplication, to an external vacuum 
solution; the azimuthal field there satisfies B — 0. At the 
inner boundary r = ro we used the same conditions as 
Tavakol et al. (2002). We note that 'reasonable' changes 
in these boundary conditions do not qualitatively change 
our conclusions (cf. Tavakol et al. 2002). 

We set ro = 0.64; with the solar CZ proper be- 
ing thought to occupy the region r > 0.7, the region 
r o < r 0.7 can be thought of as an overshoot re- 
gion/tachocline, where the major radial gradients of f2 
occur. Equations JTJ and J2J were solved using a modified 
version of the code described in Moss & Brooke (2000), 
using the the above boundary conditions, over the range 
ro < r < 1, < 6 < it. We used a mesh resolution of 
61 x 101 points, uniformly distributed in radius and lati- 
tude respectively, i.e. Ar = 0.06, A9 = 0.017T. 



2.2. Density stratification 

The density profile in the solar model computed by 
Christensen-Dalsgaard et al. (1996; henceforth referred to 
as CD96) is shown in Fig. ^ We approximated the density 
stratification in the solar convection zone by a polytropic 
density profile, writing in general 



p(r)/po 



(3) 
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Fig. 1. Density profiles obtained using a polytropic 
model, for several values of p r &tio- Here p is dimensionless, 
measured in units of the density at r — tq. The upper set 
of curves show log p, the lower show p. Shown also is the 
density profile of CD96. For p rat i 10 2 , the density pro- 
file is approximately invariant through most of convection 
zone, apart from the regions near the surface. 



where 



' l/n 
/^ratio 



'0 



l/n 
Pratio 



l 



(4) 



we take n = 1.5 below. This density profile corresponds 
to a polytrope truncated at the appropriate radius to give 
the required value of p ra tio, with this radius subsequently 
renormalized to unity. p ra tio is thus the ratio of the den- 
sities at the top (r = 1) and the bottom (r = ro) of the 
computational domain; po = p{ro). Several examples of 
such profiles, with different values of p r atio are also shown 
in Fig.n For p lat i ^ 100, these give a reasonable approx- 
imation to the density stratification in the bulk of the 
solar CZ, apart from the outermost regions near the solar 
surface 



2.3. Explicit a-quenching 

The nonlinearity in our original model is through the non- 
linearity in the Lorentz force and the subsequent presence 
of the v' term in the induction equation. In all of our 
previous studies the form of ev r was chosen to be inde- 
pendent of the magnetic field B. It did, however, include 
an implicit a-quenching in the overshoot layer in that a 
was prescribed to vanish in r < 0.7 and to be relatively 
small in the layers immediately above where \B\ was an- 
ticipated to be relatively large. Here we study the effects 
of having an additional nonlinearity in the form of an ex- 



plicit nonlinear a-quenching given (reverting temporarily 
to dimensional units) by 



a = a r (r) 



sin 



' COS f 



l+g'\B*\/B2 



(5) 



where B cq = (ATrpv 2 ) 1 / 2 is the equipartition field strength, 
v t is the turbulent velocity and g' is of order unity. (The 
trigonometric factor is chosen to agree with that used in 
our earlier work - see above.) In order to estimate how 
B^ q , appearing in the denominator, varies with depth we 
recall that to a very good approximation the luminosity 
L = Airr 2 F r oc r 2 pvt 3 is a constant in the CZ, which 
gives Vt oc / 9 _1 / 3 r~ 2 / 3 and hence pv 2 oc /7 1//3 r~ 4 / 3 . This 
suggests that the denominator in ^ will be a slowly 
varying function of the depth, especially away from the 
solar surface, and this expectation is supported by de- 
tailed mixing length models of the solar CZ. Thus we take 
B 2 q oc pv 2 — const, and so the a-quenching expression 10 
can be well approximated (now temporarily using dimen- 
sionless variables again) by 



a r (r)- 



l + g\B 2 



(G) 



where g = g' A4n por/oQ* / B 2 q . B 2 q is constant, and thus 
so is g. (In this connection it may be noted that there is 
an ongoing controversy regarding the nature and strength 
of a-quenching. This is not the place to rehearse the ar- 
guments for and against 'strong' alpha-quenching (where 
g' ^> 1); the issue is unresolved, and we use I© as a com- 
monly adopted nonlinearity.) 

We can make an estimate of what might be plausi- 
ble values for g by taking 770 = 3 x 10 11 cm 2 sec _1 n (see 
Sect. below), A = 10 4 , p = 1 gm cm" 3 , !]* = 2x 10~ 6 
sec -1 , B cq — 3500 G (the latter estimated from an en- 
velope model by Baker & Temesvary (1966) (more recent 
models do not give significantly different values). Then 
g — 3 x IOV. (As noted above, B 2 , g oc A, but B oc A~ 1/2 , 
so the results in physical variables do not depend on the 
value adopted for A.) 

The inclusion of both types of nonlinearity in the 
model allows us to study the relative importance of each 
term and in particular what happens to the amplitudes of 
the torsional oscillations as g increases. 



3. Results: effects of density stratification 

Using the above model, we studied the dynamics in the 
convection zone, and in particular the changes caused by 
the introduction of the density stratification and the ad- 
ditional nonlinearity. We concentrated on the effects of 
these changes on the torsional oscillations and in partic- 
ular on their amplitudes (i.e. the angular velocity resid- 
uals, defined as the current angular velocity minus the 
background angular velocity f2o), since these are the most 
important indicators to be compared with observations. 
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Fig. 2. The radial (r-t) contours of the angular veloc- 
ity residuals Sil as a function of time for a cut at 10 de- 
grees latitude. The model parameters are R a = —4.0, P t = 
1.0, Ru = 60000 and p ra t io = 10000. Note the penetration 
of the torsional oscillations all the way down to the bot- 
tom of the convection zone. Darker and lighter regions 
represent positive and negative deviations from the time 
averaged background rotation rate. 



The model predictions are likely to be of particular im- 
portance in deeper regions where the uncertainties in the 
inversions are largest. 

We first study the effects of density stratification in iso- 
lation, with g = 0. We fixed R u = 6 x 10 4 (corresponding 
to ?7o ~ 3 x 10 11 cm 2 sec _1 ) by requiring that the dynamo 
period is near to 22 years close to marginal stability. The 
marginal value of R a is then about —3.2. 

3.1. Comparison with previous results 

To begin with we made a detailed comparison of the re- 
sults obtained here in the presence of density stratifica- 
tion, with our previous studies which ignored this effect. 
In particular, we concentrated on one of the findings of 
these earlier studies, namely the fact that torsional os- 
cillations were found to extend all the way down to the 
bottom of the CZ. 

Our results here show that, for both near-critical, and 
moderately supercritical dynamo regimes, the oscillations 
extend all the way down to the bottom of the convec- 
tion zone, for all 1 < /O ra tio < 10 6 - This is in line with our 
previous results which shows that such behaviour is an ex- 
tremely robust feature of the models considered to density 
stratification. An example of such generic penetration is 
depicted in Fig. El This is important observationally since 
the inversions have so far not been able to resolve unam- 
biguously regimes near the bottom of the convection zone 
(Howe et al. 2000a; Antia & Basu 2000). 

For completeness we also show in Fig. [3] the butterfly 
diagram (i.e. latitude-time plot) for the toroidal compo- 
nent of the magnetic field B near the surface, as well as the 
corresponding plot for the angular velocity residuals (i.e. 
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Fig. 3. Butterfly diagram of the toroidal component of 
the magnetic field B at r — 0.99. The parameters are the 
same as those for Fig. [21 Dark and light shades correspond 
to positive and negative values of respectively. 
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Fig. 4. Angular velocity residuals at r = 0.99 with lati- 
tude and time. The parameters are the same as those for 
Figure [5] A temporal average has been subtracted to re- 
veal the migrating banded zonal flows. Darker and lighter 
regions represent positive and negative deviations from 
the time averaged background rotation rate. 



the near-surface torsional oscillations) in Fig. 21 As can be 
seen, the presence of significant density stratification does 
not alter the qualitative behaviour of the magnetic field 
and torsional oscillations. 



3.2. Amplitudes of oscillations as a function of /9 ra tio 

From an observational point of view, an important fea- 
ture of the dynamics of the CZ is the way the amplitudes 
of the torsional oscillations vary as a function of model 
ingredients and parameters, as well as depth in the CZ. 
A theoretical understanding is crucially important here, 
especially given the uncertainties in inversions of the he- 
lioseismic data in the lower parts of the CZ. 
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With this in mind, we studied how the amplitudes of 
the torsional oscillations change as a function of p Ta ±io, for 
two values of R a ; one corresponding to a near-critical and 
the other to a supercritical dynamo regime. The results 
are depicted in Figs. and For orientation, we recall 
that the corresponding amplitudes in the sun, deduced 
from helioseismological observations, are of the order of 
one nHz. 

The important trend revealed by these figures is that, 
in both the near-critical and supercritical cases, the am- 
plitudes of the torsional oscillations in the upper part of 
the CZ increase as p ra tio increases, while those at the bot- 
tom decrease. The decrease in amplitude of the torsional 
oscillations in the bottom of the convection zone, as /? ra tio 
increases (i.e. as the density profile becomes more realistic) 
is of practical importance, as it suggests that the detec- 
tion of such oscillations may be even more difficult than 
might appear at present just from the limits to the reso- 
lution of the helioseismological inversions at the bottom 
of the convection zone. The increase of the amplitudes at 
the top is also important since it improves the agreement 
between simulated and observed values. 



Also of potential interest is the way that the dis- 
tribution of the perturbation kinetic energy, defined by 
f v p{2vqv' 4- v' 2 )dV, is modified as the density stratifi- 
cation changes. In Fig. [3 we show the variation of this 
integral with radius, where we take the domain V of inte- 
gration to be the local spherical shell with thickness Ar. 
Introducing a stratification markedly reduces the pertur- 
bation kinetic energy in the near-surface regions, but the 
differences between the two stratified models illustrated 
are relatively small. However the low values of the near- 
surface density in the stratified models means that the 
amplitude of the torsional oscillations in this region are 
not correspondingly decreased - see Figs. 5 and 6. This 
relative insensitivity of the observed quantity to the strat- 
ification is reassuring. 



3.3. Amplitudes of oscillations as a function of R D 
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Fig. 5. The variation of oscillation amplitudes at r — 
0.70,0.88 and 0.95 as a function of /? r atio, for a near-onset 
dynamo regime. The parameters are R a = — 3.2, P r = 
1.0, Ru = 60000, with a — a r (r)f{6) and no a-quenching. 
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Fig. 6. The variation of oscillation amplitudes at r = 
0.70,0.88 and 0.95 as a function of p ra tio, for a supercrit- 
ical case. The parameters are R a — — 4.4, P r = 1.0, i?^ 
hOOOi). with a = a r (r)f(9) and no a-quenching. 



To begin with, we verified that for a given density profile 
(in presence of density stratification), the amplitudes of 
the oscillations increase as the dynamo number R a is in- 
creased (see also Covas et al 2001a). The results are shown 
in Fig. |S| We then studied how the amplitudes change as 
we change the density profile by increasing the value of 
Pratio- For comparison we also calculated the amplitudes 
for the density profile of CD96. The results are shown in 
Fig. As can be seen, the effect of increasing p ra tio is 
to increase the amplitudes of oscillations near the surface. 
The largest amplitudes are attained for the CD96 density 
distribution. 



3.4. Amplitudes of oscillations as a function of depth 

and /3 ra tio 

We also studied the variations of the amplitudes of oscil- 
lations as a function of depth for a number of values of 
Pratio for both a near-onset regime (with R a = —3.2) as 
well as for more supercritical parameters (R a — —4.4), 
and the results are shown in Figs. 1101 andllll respectively. 
As can be seen, despite the latter having larger ampli- 
tudes, the two cases have important similarities. In each 
the amplitudes near the surface layers increase with in- 
creasing Pratio, whereas they decrease near the bottom. 
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Fig. 7. The distribution of the kinetic energy (in arbitrary 
units) in spherical shells of thickness the radial mesh size 
Ar, as a function of radius and p ra tio, for both near-onset 
and supercritical dynamo regimes. The parameters are as 
in Figs. and respectively. 



0.4 




0.0 — 1 — ' — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — ' — 1 — 1 — 1 
-3.2 -3.4 -3.6 -3.8 -4.0 

R„ 

Fig. 9. The variation of oscillation amplitudes as a func- 
tion of R a for different values of p r atio at r — 0.95. Shown 
also are the variations obtained using the solar model den- 
sity profile of CD96. The parameters are R u = 60000, 
P r = 1.0 with a = a r (r)f(9) and no a-quenching. 
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Fig. 8. The variation of oscillation amplitudes as a func- 
tion of R a . Parameter values are P r = 1.0 and i? w = 
60000, pratio = 10000, with a = a r (r)f(0) and no a- 
quenching. 

4. Results: comparative study of the nonlinearities 

In order to make a comparative study of the effects of the 
two competing forms of nonlinearities considered here, we 
begin (see also below) by recalling that, in the absence of 
any alpha-quenching, increasing the degree of nonlinearity 
due to the Lorentz force (i.e. increasing R a ) results in an 
increase in the amplitudes of the torsional oscillations. The 
question is what is the effect of increasing the strength of 
the additional nonlinearity due to a-quenching on these 
amplitudes? 

We studied this question by making a detailed study 
of how the amplitudes of torsional oscillations change as 
a function of g, the coefficient of the a-quenching term 



0.12 




0.00 1 

0.6 0.7 0.8 0.9 



r 

Fig. 10. The variation of oscillation amplitudes as a func- 
tion of depth for a number of values of /? ra tio, for a near- 
onset regime with R a = -3.2, R u = 60000, P r = 1.0, 
a = a r (r)f(9) and no a-quenching. 

(Eq. (JSJ). The results are presented in Fig. ^| A num- 
ber of trends are apparent. The amplitudes of the tor- 
sional oscillations decrease at all depths with increasing 
a-quenching coefficent g. However, the oscillations at the 
top of the CZ have the highest amplitudes throughout. 

To study the influence of increasing the a-quenching 
coefficient, we recall that in the case of pure alpha- 
quenching, 

< B 2 >oc (7) 
9 

where the angled parentheses denote some sort of global 
average. We checked that this was exactly satisfied for 
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Fig. 11. The variation of oscillation amplitudes as a func- 
tion of depth for a number of values of p ra tio, for a super- 
critical case with R a = -4.4, = 60000, P r = 1.0, 
a = a r (r)f(9) and no a-quenching. 

our code in the quasi-kinematic approximation when the 
Navier-Stokes equation was switched off. Interestingly, we 
also found similar behaviour when the additional nonlin- 
ear feedback via the Navier-Stokes equation was present. 
This is shown in Fig. which displays the behaviour of 
the global magnetic energy as a function of g. 

The amplitudes of the torsional oscillations increase 
somewhat with \R a \ but in the range accessible to our 
code (\R a | Si 15) this increase is insufficient to compensate 
for the decrease caused by the increase in the a-quenching 
coefficient g (given by (7)). 

We can add here that in our previous, very prelimi- 
nary, study of the effects of alpha-quenching on the tor- 
sional oscillations (Covas et al 2001a), we were only able 
to find oscillatory solutions for very small values of the 
alpha-quenching coefficient (g < 0.1). What has enabled 
us to find the existence of oscillatory solutions at such 
large values of g here is our recent discovery of the ex- 
istence of multiple attractors in the dynamo system, in- 
cluding both oscillatory and steady solutions, that coexist 
to g ~ 1000, i.e. in the physically plausible range as esti- 
mated in Sect. 2.3. 

5. Discussion 

We have made a detailed study of the effects of including 
density stratification in the model, as well as an additional 
form of nonlinearity due to a-quenching on the behaviour 
of solar torsional oscillations. We feel that the previous 
omission of any discussion of these effects reduced the de- 
tailed applicability of our results to the observational data. 

To study the effects of density stratification, first of all 
we made a comparison with our previous results and found 
both important similarities and differences. Concerning 
similarities, we found that consistent with our previous 
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Fig. 12. The variation of amplitudes of torsional oscilla- 
tions as a function of changes in the a-quenching coeffi- 
cient 3 for a near-onset regime. Parameters are R a = —2.5, 
P T = 1.0, Ru; = 60000 with a r = 1 (i.e. no explicit radial 
dependence of a) and the density p given by the solar 
model of CD96. 
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Fig. 13. The variation of the magnetic energy as a func- 
tion of the a-quenching coefficient g. Shown arc both 
the kinematic (without the Navier-Stokes equation) and 
the nonlinear (with the Navier-Stokes equation) cases. 
Parameters are R a = —2.5, P r = 1.0 and R w = 60000, 
with a r = 1 and the density p given by the solar model of 
CD96. 



findings for near-critical and moderately supercritical dy- 
namo regimes, the torsional oscillations extend all the way 
down to the bottom of the convection zone. This is true in 
presence of density stratification with all values of p r atio in 
the range [1, 10 6 ]. Concerning differences, we found that 
in both the near-critical and supercritical cases, the am- 
plitudes of the torsional oscillations in the upper part of 
the CZ increase with Pratio while those at the bottom de- 
crease. Given that the uncertainties in the inversions are 
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expected to be greatest at the bottom of the CZ, the lower- 
ing of amplitudes in these regions with increasing density 
stratification has potential importance for the analysis of 
helioseismological data, and for the understanding of the 
dynamo regimes near the bottom of the convection zone. 

When two forms of nonlinearity are simultaneously 
present, we find that torsional oscillations persist in pres- 
ence of both nonlinearities, for a large range of values of 
the a-quenching coefficient g up to g > 1000. The am- 
plitudes of these oscillations decrease however, as they 
shadow the square of the magnetic field strength, which 
in turn follows the scaling with g given by (JJJ. 

We also note that for small g, the combined effect of 
increasing /0 rat i o and g is to leave the amplitude of the near- 
surface torsional oscillations little changed. However, for 
larger g (> 10), the effect of a-quenching dominates, as 
the amplitudes of the torsional oscillations scale as 1/g, 

Equation J5J suggests, and our models confirm, that 
variations in v' closely shadow those in the ^-component 
of the Lorentz force. It seems a priori difficult to imag- 
ine a driver for torsional oscillations that does not in- 
volve the Lorentz force, either at a macro- or micro-level 
(e.g. the 'nonlinear A-effect', Kitchatinov & Pipin 1998; 
Kitchatinov et al. 1999), whatever the prime nonlinearity 
in the dynamo. Thus the requirement that the Lorentz 
force shadows the torsional oscillations is plausibly an ad- 
ditional constraint even on dynamo models that do not 
explicitly include dynamical effects (e.g. conventional a- 
quenched models. Brandenburg & Tuominen 1988 appear 
to have been the first to investigate semi-quantitatively 
such an idea.) However it is perhaps a little puzzling 
that values of g as large as suggested by the estimates 
of Sect. 12.31 reduce the amplitudes of the torsional oscil- 
lations so substantially, removing the previous agreement 
with observations. 
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